function [feature] = caculate_feature(filename)
    fs = 100;
    fftpoint = 256;
    
    window = 100;
    overlap = 0;
    thd_alpha = 25;
    thd_beta = 25;
    thd_theta = 25;
    thd_delta = 50;
    data = load(filename, '-ascii');
    c = xcov (data);
    L = length(c)
    
    Y = fft(c,fftpoint)/L;

    L = fftpoint;
    al_L = L /2
    POWER = 2*abs(Y(1:L/2+1));
    f = fs/2*linspace(0,1,L/2+1);
    %frequency
    segment_point0 = 0.6 * fftpoint/fs;
    segment_point1 = 4 * fftpoint/fs;
    segment_point2 = 8 * fftpoint/fs;
    segment_point3 = 13 * fftpoint/fs;

    alpha = POWER(segment_point2:segment_point3); 
    beta = POWER(segment_point3:L/2+1);
    theta = POWER(segment_point1:segment_point2); 
    delta =  POWER(segment_point0:segment_point1); 
    
    totalpower_alpha = sum(alpha);
    totalpower_beta = sum(beta);
    totalpower_theta = sum(theta);
    totalpower_delta = sum(delta);
    totalpower = totalpower_alpha + totalpower_beta + totalpower_theta + totalpower_delta;

    averagepower_alpha = totalpower_alpha/length(alpha);
    averagepower_beta = totalpower_beta/length(beta);
    averagepower_theta = totalpower_theta/length(theta);
    averagepower_delta = totalpower_delta/length(delta);
    averagepower = averagepower_alpha + averagepower_beta + averagepower_theta + averagepower_delta;

    per_aver_alpha = averagepower_alpha / averagepower;
    per_aver_beta = averagepower_beta / averagepower;
    per_aver_theta = averagepower_theta / averagepower;
    per_aver_delta = averagepower_delta / averagepower;
    %rate
    delta_alpha = averagepower_delta / averagepower_alpha;
    delta_beta = averagepower_delta / averagepower_beta;
    delta_theta = averagepower_delta / averagepower_theta;
    theta_alpha = averagepower_theta / averagepower_alpha;
    theta_beta = averagepower_theta / averagepower_beta;
    alpha_beta = averagepower_alpha / averagepower_beta;
    deltatheta_alphabeta = (averagepower_theta + averagepower_delta) / (averagepower_beta + averagepower_alpha);
    
    %time
    data_mean = mean(data);
    data_variance = var(data);
    data_rms = rms(data);
    data_skewness = skewness(data);
    data_kurtosis = kurtosis(data);
    
    % time frequency
    s = spectrogram(data, window, overlap, fftpoint,fs)
    delta_point = segment_point0 / 2;
    theta_point = segment_point1 / 2;
    alpha_point = segment_point2 / 2 ;
    beta_point = segment_point3 / 2;
    delta_time = sum((sum(abs(s(delta_point : theta_point,:))> thd_delta))');
    theta_time = sum((sum(abs(s(theta_point : alpha_point,:))> thd_theta))');
    alpha_time = sum((sum(abs(s(alpha_point : beta_point,:))> thd_alpha))');
    beta_time = sum((sum(abs(s(beta_point:al_L,:))> thd_beta))');
    % localmin_max analysis
    delta_band = abs(sum(s(delta_point : theta_point,:)))
    theta_band = abs(sum(s(theta_point : alpha_point,:)));
    alpha_band = abs(sum(s(alpha_point : beta_point,:)));
    beta_band = abs(sum(s(beta_point :al_L,:)));
    
    local_min_max_delta = local_min_max(delta_band);
    local_min_max_theta = local_min_max(theta_band);
    local_min_max_alpha = local_min_max(alpha_band);
    local_min_max_beta = local_min_max(beta_band);
    
    % entropy
    %normalized
    total_data = sum(POWER);
    normalized_power = POWER / total_data;
    entropy_power = - 1 / log(length(POWER)) * sum(normalized_power .* log(normalized_power));
    
    feature = [averagepower, per_aver_alpha, per_aver_beta, per_aver_theta, per_aver_delta, ... 
        averagepower_alpha, averagepower_beta, averagepower_theta, averagepower_delta, ...
        delta_alpha, delta_beta, delta_theta, theta_alpha, theta_beta, alpha_beta, deltatheta_alphabeta, ...
       data_mean, data_rms, data_variance, data_skewness, data_kurtosis, ...
       alpha_time,beta_time, theta_time, delta_time, ...
       local_min_max_alpha, local_min_max_beta, local_min_max_theta, local_min_max_delta];
    
   
    
   
    
    
end